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ABSTRACT 

We develop a numerical hydrodynamics code using a pseudo-Newtonian formulation 
that uses the weak field approximation for the geometry, and a generalized source 
term for the Poisson equation that takes into account relativistic effects. The code 
was designed to treat moderately relativistic systems such as rapidly rotating neutron 
stars. The hydrodynamic equations are solved using a finite volume method with High 
Resolution Shock Capturing (HRSC) techniques. We implement several different slope 
limiters for second order reconstruction schemes and also investigate higher order re- 
constructions such as PPM, ENO and WENO. We use the method of lines (MoL) to 
convert the mixed spatial-time partial differential equations into ordinary differential 
equations (ODEs) that depend only on time. These ODEs are solved using second 
and third order Runge-Kutta methods. The Poisson equation for the gravitational 
potential is solved with a multigrid method, and to simplify the boundary condition, 
we use compactified coordinates which map spatial infinity to a finite computational 
coordinate using a tangent function. In order to confirm the validity of our code, we 
carry out four different tests including one and two dimensional shock tube tests, sta- 
tionary star tests of both non-rotating and rotating models and radial oscillation mode 
tests for spherical stars. In the shock tube tests, the code shows good agreement with 
analytic solutions which include shocks, rarefaction waves and contact discontinuities. 
The code is found to be stable and accurate: for example, when solving a stationary 
stellar model the fractional changes in the maximum density, total mass, and total an- 
gular momentum per dynamical time are found to be 3 x 10 -6 , 5 x 10 -7 and 2 x 10 -6 , 
respectively. We also find that the frequencies of the radial modes obtained by the 
numerical simulation of the steady state star agree very well with those obtained by 
linear analysis. 

Key words: relativistic processes - gravitation - hydrodynamics hydrodynamics- 
methods: numerical 



1 INTRODUCTION 

It is necessary to take into account both special- and general 
relativistic effects in the studies of the dynamics of com- 
pact astrophysical object such as neutron stars and black 
holes. Some pulsars produce pulses of up to 1 KHz, corre- 
sponding to rotation speeds at the surface of around 0.2c. 
Their typical sizes and masses are known to be around 
10km and 1.4 ~ 2M0, respectively, giving compactness, 
GM/Rc 2 = 0.2 ~ 0.3. Therefore, a Newtonian approach 
cannot properly describe neutron stars, even for the non- 
rotating case. 
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In general relativity, the dynamics of gravity (or space- 
time) can be studied by solving the Einstein equations. 
The equations of motion for the matter are given, in part, 
by the conservation law of the energy-momentum tensor 
which itself sources the gravitational field. Computational 
approaches for solving general relativistic field equations 
constitute the field of numerical relativity. 

Over the past few decades, many general relativis- 
tic hydrody n amic codes have been developed, starting 
with IWilsonl (|l972t) who proposed a 3+1 Eulerian formu- 
lation (see also IWilson fc Mathewl 120031 '). Although Wil- 
son's numerical approach was widely used to study prob- 
lems such as core collapse and accretion disks, it pro- 
duced large errors when fluid flows became ultra-rela tivistic 
ijCentrella fc Wilson|[l983 : iNorman fc Winklerlll98rj) . In or- 
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der to avoid these excessi v e erro rs, a new formulation was 
proposed by iMarti et alj |l99l[ ). This formulation makes 
it possible to use existing numerical techniques based on 
characteristic approaches for Newtonian hydrodynamics. In 
particular, these include High Resolution Shock Capturing 
(hereafter HRSC) methods that reduce the order of accuracy 
near shocks, but minimize the amount of numerical dissipa- 
tion. This dissipation is very unnatural and can result in 
non-physical effects in the numerical results. Marti's formu- 
lation was exte nded to the gene ral relativistic case by the 
Valencia group l|Font et al.ll200oh . and this last work forms 
the basis for most recent general relativistic hydrodynami- 
cal codes. Recent reviews of the formulation and num erical 
metho ds can be found in IMarti fc Mulled <|2003h and iFontJ 
<l200Sh . 

However, when working in multiple spatial-dimensions, 
it still requires a lot of computational resources to treat fluid 
dynamics in concert with the evolution of the general rela- 
tivistic gravitational field. In addition, numerical relativity 
simulations have frequently encountered instabilities which 
are often associated with violations of the Hamiltonian and 
momentum constraints. (However, with the development of 
new formulations which cast the Einstein equations in ap- 
propriate hyperbolic forms, as well as the use of constraint- 
damping tech niques, significant progre ss has been made on 
this front: see ISarbach fc Tiglio 1120121 for a very recent re- 
view of this subject). For these reasons, simulations using 
Newtonian gravity are still used even though they are not 
applicable to very compact objects. 

The aims of this paper are 1) to introduce a new 
formulation whi ch applies a pseudo-Newtonian approach 
ijKim et al.ll2009h to the study of moderately relativistic ob- 
jects and 2) to describe a numerical implementation of this 
method. In our pseudo-Newto nian approach, which was in- 
troduced bv iKim et al.l ((2009) for steady state models, the 
gravitational field is treated by a weak field approximation, 
but special relativistic effects are correctly taken into ac- 
count. Specifically, the Newtonian gravitational potential 
that appears in the weak field metric satisfies a Poisson 
equation, but the mass density that appears as a source 
term for that equation is modified to include relativistic ef- 
fects. Of course this m ethod cann o t be applied to highly 
relativistic systems, but IKim et all (|2009l ) showed that the 
pseudo-Newtonian formulation is valid for the modeling of 
mildly compact objects, such as rotating neutron stars hav- 
ing su rface rotation ve locity up to ~ 0.2c and compactness 
~ 0.2 (|Kim et all2008l ). In this paper, we extend the pseudo- 
Newtonian approach to hydrodynamical systems where the 
flows can be ultra-relativistic and gravity can be moderately 
strong. 

The remainder of the paper is structured as follows: In 
section [2] we present the formulation and governing equa- 
tions for our system, while the numerical techniques em- 
ployed in our study are given in section[3] We discuss various 
numerical tests of our code's treatment of hydrodynamics for 
the case of shock tubes in [4] and for stationary stars in [5] A 
test which compares radial pulsation mode frequencies for 
polytropic stars determined through dynamical evolution to 
those computed in linear theory is detailed in section [6] We 
conclude with a summary and discussion in section [7] 

Throughout this paper we use units in which c = G = 



Mq = 1: these correspond a to unit time = 4.92 x 10 3 ms, 
unit length = 1.47km and unit mass = 1.99 x 10 33 g. 



2 FORMULATION 

Our pseudo-Newtonian method was first discussed in the 
steady state context by IKim et all l|2009h . We assume the 
weak field metric, 

ds 2 = g^dx^dx" 

= -(1 + 2<S>)dt 2 + (1 + 2§y 1 8 i jdx l dx j , ' ' ' 

where g M „ is the spacetime metric and 3> is the Newto- 
nian gravitational potential. With this metric, we neglect 
all higher order effects such as frame dragging and describe 
gravity using only a single gravitational potential, just as in 
the Newtonian case. The gravitational potential satisfies a 
Poisson equation with the active mass density, p a ctivc, pro- 
viding the source: 

V 2 $ = 47Tpactive. (2) 

The active mass density is computed from the relativistic 
definition of the energy. For a perfect fluid, the energy mo- 
mentum tensor can be expressed as 

T" u = p hu^u v + Pg^ ', (3) 

where the specific enthalpy is defined by 

h = l + e+— . (4) 
Po 

The active mass density is then given by 

Pactivc = T — 2T ° = Tt — r ° = pohl^- + 2P. (5) 

1 — tr 

In Eqs. 0, Q and ([5]), po is the rest mass density which is 
proportional to the number density of baryons of the fluid, 
P is the pressure, it M is the four velocity of a fluid element 
with respect to an Eulerian observer, e is the specific internal 
energy, T is the trace of the energy-momentum tensor (T = 
g^vT^), and v is the three dimensional fluid velocity. Unlike 
the Newtonian case, Pactive includes all sources of energy. 

The equations governing the motion of the fluid mat- 
ter can be derived from the conservations laws for the 
energy-momentum tensor and the fluid's matter current, i.e., 
V M T M " = q anc j y u jM = q In ^ he ADM decomposition of 

spacetime jArnowitt et alj|l962f ). the metric (g M „) can be 
expressed in the following form by considering the foliation 
of spacetime using 3-dimensional hypersurfaces defined by 
t — const.: 

ds 2 = -adt 2 + 7^ (dx i + p l dtj (dx j + f3 j dt\ , (6) 

Here jij is the spatial metric, defined on each hypersurface, 
while a and /3 l are known as the lapse, and shift vector, 
respectively, and encode the 4-fold coordinate freedom of 
general relativity. 

Flux-conservative formulations of hydrodynamics have 
been applied very successfully in computational fluid dy- 
namics. To cast the fluid equations in flux-conservative form 
we first define so-called conservative variables (q) in terms 
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of the original hydrodynamic variables (so-called primitive 
variables, w), 
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where W — 1/ y/l — jijV i vi . With these definitions, and 
with th e metric (E 
tion as (|Font et al 



with th e metric (Eg . (O) , we can then write the Euler equa- 
l.l l2OO0l l 
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dt dx i 
where the fluxes f and the sources E are given by 
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Here y^y and \f—g are the determinants of jij and g M „, re- 
spectively, and are related by y/—g = otyFy. It is well known 
that for a perfect fluid, the system of equations derived from 
the conservation laws is not closed: the number of dynamical 
equations is always less than the number of unknowns. 

As is also well known, the equation of state (hereafter 
EOS) for the fluid provides an additional equation, but in 
the general case it also introduces other unknowns. In order 
to completely close the hydrodynamical equations, an en- 
ergy balance equation is often used. However, under certain 
circumstances, we can adopt rather simple EOSs that do 
not introduce any further variables: adiabatic and isother- 
mal EOSs provide specific examples. 

Realistic EOSs are usually determined by theoreti- 
cal calculations and experimental measurements. However, 
there are physical regimes where our understanding of the 
nature of the matter is quite incomplete. Specifically, in the 
case where the matter density is significantly above nucleon 
density, there remain large uncertainties in the correct EOS. 
Thus, for example, the EOS at the core of neutron stars is 
still not very well understood. Here we ignore these difficul- 
ties, and for the purpose of testing our code, use two types 
of very simple EOS. The first is the ideal gas EOS which 
can be written in the following form: 



P=(r-l)p e, 



(10) 



and corresponds to the isothermal EOS. We use this EOS 
in the shock tube tests described in (section S]). The second 
EOS results from the isentropic assumption, whereby Eq. 
pop becomes the polytropic EOS: 



P = Kp 



1+- 



(11) 



Here K and N are the polytropic constant and index respec- 
tively. The polytropic EOS of state is the generalized form 
of the adiabatic one; a fluid which is governed by it does 
not generate entropy, and shock formation is thus generi- 
cally prohibited. We use this EOS in the pulsation mode 
test (section [6] and Appendix |A"|) . 



Using the above formulation, we are now ready to de- 
scribe in detail the pseudo-Newtonian hydrodynamical equa- 
tions used in our code. We limit our study here to axisym- 
metric systems, and adopt cylindrical coordinates (R, Z, <fi) 
such that 

ds 2 = -(1 + 2$)di 2 + i 2$ (dR 2 + dZ 2 + R 2 d<f> 2 ) . (12) 

The lapse function and shift vector are thus given by a = 
VI + 2$ and /3 l = 0. In addition, we enforce the equatorial 
symmetry at 2 = since the phenomena involving I = odd 
modes are not dominant in most cases of neutron star dy- 
namics, where I is from the spherical harmonics. In this coor- 
dinate system, the conservative and primitive variables are 
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The final form of the hydrodynamical equations then be- 
comes 



Of 
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Using Eq. (JT2J we have ^7 = ^(1 + 2$)" 3/2 , and ^fg = 
R (1 + 2"I>) _1 . In obtaining the expressions in Eq. (|15|) we 
have used the assumption of slow changes of the poten- 
tial rel ative to the gr a dients (^- <g ^ or §§)• Re- 
cently, iNagakura et al.l (|201ll ) used a similar method in 
the context of jet propagation in a uniform medium, but 
adopted a slightly different linear momentum equation than 
ours. (See Eq. ( | 15t and compare with Eqs. (2) and (3) in 
INagakura et al-lHoill ). 

Finally, the gravitational Poisson equation in our coor- 
dinate system is 
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(16) 



Note that the second component of E contains terms 
which, individually, become singular on the axis of symmetry 
(R — 0). In addition, there are other terms in the equations 
of motion that need to be treated carefully as 7? — > 0. This is 
done by demanding regularity at the axis, and by considering 
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the parity of each function, with respect to R, in that limit. 
In particular, po, v z , if , P, D, Sz, S$, and r are all even 
functions of R as R — » 0, while v fl and Sr are odd. Taking 
this into account, Eq. (|14[) and £ in Eq. (|15[) become 



and 



E = 



1+2$ BZ 
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(18) 



(1 + 2<S>) _1 . The 



in Eq. (|17[) becomes 2 instead of 1, 



where V7 7 = (l + 2$)" 3/2 and 

coefficient of 8 ^ v ^^ — - 
while the other variables, such as q, f R and / , are un- 
changed from Eq. (fl5)l . 

Finally, using L'Hopital's theorem at R = 0, the singu- 
lar term R~ x d<& / dR in the Poisson equation (|16[1 is replaced 
by d 2 <£>/dR 2 . 



3 NUMERICAL METHODS 

In this section, we describe our numerical methods for solv- 
ing the coupled hydrodynamical and Poisson equations. We 
mainly use the finite volume methods for the hydrodynam- 
ical equations and the a finite difference approach for the 
Poisson equation. In the finite volume method, each grid cell 
represents volume averaged hydrodynamic quantities i.e., 
q = -^y J qdV. After applying the finite volume method our 
hydro equations can be reduced to Riemann problems which 
consider the time evolution of initial conditions given by two 
distinct states that adjoin at some interface (so that there 
are, in general, discontinuities across one or more physical 
quantities at the interface). A very important property of 
the finite volume method is that it maintains the local con- 
servation properties of the flow in the computational grid. 

In the dynamics of compressible fluids, we inevitably 
encounter discontinuous behaviors such as shocks, rarefac- 
tions or contact discontinuities. To treat such discontinuities 
without introducing numerical instabilities or spurious os- 
cillations, we use High Resolution Shock Capturing (HRSC) 
techniques that generically reduce the order of accuracy of 
the numerical scheme near discontinuities or when one or 
more of the fluid variables are at a local maximum. A key 
ingredient to the success of the HRSC methods is the calcu- 
lation of fluxes through cell boundaries. To compute these 
fluxes we need approximate values for the primitive vari- 
ables at the cell boundaries. We ha ve implemented second 
order slope limiters such as minmod (|van Leerl ll979). mono- 
tonized ce ntral differe nce (MC hereafter. Ivan Leerll 19771 ) and 
superbee (|Roel 1 1985ft. as well as a third order slope lim- 
iter proposed by Shi bata (2003) and which is based on the 
minmod function (3minmod hereafter). Other reconstruc- 
tion methods such as the third order Piecewise Par abolic 
Method (PPM hereafter, IColella fe Woodward! I1984T ), Es- 



sentially Non-Oscillatory meth od fENO lHarten et alj 1987) 
and W eighted ENO (WENO, iLiu et al.|[l*994l ; Ijiang fc Shu] 
1996), which has an arbitrary order of accuracy, were also 
implemented. 



In the implementation of HRSC schemes it is not ef- 
ficient to exactly solve the Riemann problems which arise 
since an excessively large amount of computational resources 
per cell are then needed to calculate the fluxes. Thus, an 
approximate calculation of fluxes is perfo rmed. We imple- 
mented J;heJ)o21owirig_tlrr^ Mar- 
quina (iDonat fc Marquhial Il996l: [bonatl Il998ft, and H LLE 



jHarten et al.lll983l ; lEinfeldtlll988l ; lEinfeldt et ailll99lft ap- 
proximations. The Roe approximation is based on the 
Rankine-Hugoniot jump condition and Marquina's approach 
generalizes Roe's scheme. The HLLE algorithm comes from 
a very simple two wave approximation and produces the 
most dissipative and stable results. We have mainly used 
HLLE for reducing computational cost but found that our 
results were not significantly influenced by the type of flux 
approximation used. 

In order to solve the Poisson equation (which is el- 
liptic) for the gravitational potential we use the multigrid 
method, which can quickly reduce low frequency error com- 
ponents in the solution by adopting hierarchical grid levels 
|Brandd [l977ft . One of the difficulties we often encounter 
with the Poisson equation is in the proper implementation 
of the boundary conditions. For example, one of the natural 
boundary conditions is $ = at 00, but in the coordinates 
adopted in the previous section the computational domain 
cannot reach spatial infinity. We thus now refer to our pre- 
vious coordinates as (r, z) and introduce new coordinates 
(R, Z) which compactify the spatial domain, mapping the 
infinities in each spatial direction to finite coordinate values. 
Specifically, we choose the same type of compactification for 
both r and z coordinates, namely a tangent function, but 
allow a certain portion of the domain to remain "'uncom- 
pactified" : 



R if R < r 

r + n tan (^ ra ) i£R>r 

Z if Z < z 

z + 21 tan (-^t 2 -) if Z > z 



(19) 
(20) 



Here, the four parameters zq, Zi, ro and n control the com- 
pactification, and we chose this specific form for the co- 
ordinate transformation since it guarantees that the com- 
pactified coordinates smoothly transition to the original 
ones near the origin. We note that we solve the hydrody- 
namical and gravitational equations on separate spatial do- 
mains: [0 : ro,0 : Zo] for the hydrodynamic calculations and 
[0 : ro + -ri, : zo + —zi] for the computation of the gravi- 
tational potential. There ranges correspond to [0 : ro, : Zq] 
and [0 : co, : 00], respectively, in the original cylindrical co- 
ordinates (r, z). In the compactified coordinates, the Poisson 
equation is written as 



rf(R) dR 
1 d 



+ - 



r d$(R, Z) 
f(R) dR 
1 d®(R, Z) 



g(Z) dZ \g{Z) dZ 



— 47rp ac ti„ 



(21) 
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where f(R) and g(Z) are given by 



/(*) = 



g(z) = 



dr 1 


r 1 


dR ~ | 


sec^ 
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f 1 , 




sec^ 



.2 f Z-zq 
*1 



if i? < r 
if R > r , 

if Z < z 
\i Z > z a 



Table 1. Initial values of physical quantities for Shock tube tests 
(Riemann problem). 



(22) 
(23) 



As just noted, the domain for the hydrodynamical calcula- 
tion is finite, i.e. we do not solve the hydrodynamical equa- 
tions on the full compactified domain, and we thus must be 
careful to choose values of ro and zo large enough so that 
there is no outfiux of matter through the r = ro and/or 
z — zo boundaries. In our code we set ro = zo = We 
where r\ is a free parameter and r e is the equatorial ra- 
dius of the rotating star as obtained from the procedure we 
use to calculate the initial stellar model. For the pulsation 
mode test described in section [(J] a typical choice is r\ = 2. 
This means that the hydrodynamical computational domain 
extends twice the distance of the stellar radius in both the 
R and Z directions: this choice is found to be sufficient for 
our study. The values of ri and z\ are automatically de- 
termined by requiring the multigrid domain to be 2 times 
larger than the size of hydrodynamic domain in compactified 
coordinates, i.e., ri = -|ro and z\ = \zo. 

In our multigrid algorithm, we use line relaxation for 
our basic smoother, whereby all grid point values given 
by R — const, or Z — const, are updated simultaneously 
(constant- J? and constant- Z sweeps are alternated). We can- 
not use point-wise relaxation since, as is well known, such 
a technique is not a good smoother when there is signif- 
icant anisotropy in the coefficients of the second deriva- 
tive terms in the elliptic operator being treated. This is 
the case in our compactified coordinate system, particu- 
larly near the domain boundaries. We have used second- and 
fourth-order finite-difference approximations to the Poisson 
equation, and these lead to tri- and pentadiagonal linear sys- 
tems, respectively, that must be solved to implement the line 
relaxations. We use the routines DGTSV (tridiagonal) and 
DGBSV (banded/pentadiagonal) routines from LAPACK to 
perform these solutions. 

In order to integrate the discretized hydrodynamical 
equations, we use the method of lines (MOL), transform- 
ing our partial differential equations in time and space to 
ordinary differential equations (ODEs) with respect to the 
time. To solve these ODEs, we then employ second and third 
order Runge-Kutta methods, which are known to have the 
Total Variation Diminishing (TVD) property. 



4 SHOCK TUBE TESTS 

In order to verify the accuracy and the convergence of our 
numerical code, we first carry out rigorous tests using initial 
configurations having analytic solutions. In this section, we 
present the results of such tests for the case where there 
is no self-gravity (i.e. pure hydrodynamics). Another test of 
the entire code — including our treatment of the gravitational 
field — is described in the next section. 

Shock tube tests are Riemann problems where the ini- 
tial configuration of the fluid is given by two states having, 
in general, different densities, pressures and velocities, on 
the left and right halves of the tube. Three possible distinct 
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features emerge from the subsequent evolution: a shock, a 
rarefaction fan, and a contact discontinuity. We carried out 
ID and 2D numerical si mulations with 3 different pa rameter 
sets previously used by IZhang fc MacFadvenl lj200rj ). These 
parameters are listed in Tabled where superscripts R and 
L represent the fluid states in the right and left halves, re- 
spectively, of the tube. 



4.1 ID test in Cartesian Coordinates 

We first we present the results of our ID tests. In these tests, 
we carefully examine how the distinct features predicted by 
the analytic solutions are reproduced by different methods, 
and measure the accuracy and convergence rate of the vari- 
ous solutions obtained. 

In problem 1, the initial discontinuity gives 3 different 
types of solutions (shock, rarefaction, contact discontinu- 
ity). Figure 1 shows the results at t — 0.4 obtained using 
four different methods of reconstruction: minmod (top left), 
MC (top right), 3minmod (bottom left) and PPM (bottom 
right). We observe that the minmod method is quite dissipa- 
tive, yielding rather smooth solutions that cannot accurately 
describe the Shockwave. We also find that at low resolution 
the height of the shock is not well reproduced if we use the 
minmod method. MC and 3minmod give almost similar re- 
sults, while PPM shows the best behaviour near the shock. 

The second test problem (Problem 2) is the so-called 
blast wave test which produces a very sharp and thin shell 
in density between the shock and contact discontinuity. Gen- 
erally, numerical codes are not able to perfectly resolve this 
very thin shell because it can span only a few grid cells, 
even in very high resolution calculations. Nonetheless, this 
test provides insight as to how well a code can handle such a 
feature. As can see in Figure |4~T1 PPM again gives the best 
results, although it still shows large errors at the shock. 

The third problem generates a strong reverse shock 
but numerical solution has oscillatory features near the 
shock front. Generally speaking, the oscillation can be easily 
damped out if the numerical scheme is significantly dissipa- 
tive. Numerical dissipation also tends to weaken the sharp- 
ness of the discontinuity. In Figure |4~T1 one can see that the 
minmod methods, which, as already noted, is the most dis- 
sipative of the techniques we use, gives relatively small am- 
plitude oscillations, except near the discontinuity. The more 
non-dissipative methods describes the shock features well, 
but produce rather large amplitude oscillatory behavior. 

To quantify the deviation of our numerical results from 
the analytic solutions, we use the L\ norm of the errors, de- 
fined by L\ — X)i=i Axi\Qi — <l{ x i)\i where q(xi) is the value 
of the analytic solution at point Xi. We summarize the re- 
sults in Table 14. ll The convergence rate (log 2 [L\ h /L\] ) in 
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Figure 1. One dimensional shock tube test of problem 1 at t = 0.4 with different reconstruction methods: minmod(top-left), MC (top 
right), 3minmod (bottom left), and PPM (bottom right). The initial discontinuity is at x = 0.5. We use 5f2 uniform grid points. The 
numerical results are shown in 3 different colours: rest mass density (pink), pressure (red) and velocity (blue). The solid lines show the 
analytic solutions. 
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Figure 2. Same as Fig. 1 for problem 2. 
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Figure 3. Same as Figs. 1 and 2 for problem 3. 



the table should be close to 1, which corresponds to the 
1 st order nature of the HRSC scheme near the shock where 
the most of the Li-norm error occurs. However, it can devi- 
ate from that value due to the oscillatory features near the 
shock. 

Figure 0] shows the L\ norms and convergence rates for 
each problem when the grid resolution is N = 512. Although 
no single method stands out in our ID shock tube tests, 
we conclude from from the values of the L\ error norms 
and convergence rates that PPM gives the most promising 
results. 




o.oi 



MC 



4.2 2D test in Cylindrical Coordinates 

Since the cylindrical coordinate system we have adopted is 
curvilinear, 1-dimensional shock tube tests are not sufficient 
for assessing our code's accuracy and convergence. In Carte- 
sian coordinates, fluxes between cells which have the same 
state cancel out. For example, if we carry out the shock 
tube test in the x-direction, then the fluxes in the y and z 
directions are identical in every grid cell, meaning that the 
net flux is 0. Therefore, ID shock tube tests performed with 
codes that use 2- or 3D Cartesian coordinates produce ex- 
actly the same results as a ID code. However, in cylindrical 
coordinates, fluxes do not cancel in this way, but rather are 
balanced by source terms. This difference may give addi- 
tional non-physical effects, especially near discontinuities. 

Therefore, we carried out the first of the shock tube 
tests listed in Table [4] in cylindrical coordinates, where we 
placed the discontinuity on the Z—0 plane. Figure 5 shows 
the solution resulting solution on the Z-axis. If we use the 
minmod method, the 2D results are similar to the 1 dimen- 
sional ones. In addition, although PPM produces better re- 



O 0.6 
0.1 



1 


1 1 1 

Problem 1 ■ 






Problem 2 1 






Problem 3 1 






i i 





minmod MC 3-minmod PPM 



Figure 4. The L\ norm (top) and convergence rate (bottom) 
when the number of grid points is 512 with different reconstruc- 
tion methods. Three different shock tube problems are shown with 
different colors (problem 1: blue, problem 2: red and problem 3: 
sky blue). 
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Table 2. The L\ norm of the error and its convergence rate for each of the test problems using different resolutions and different 
reconstruction schemes. 
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Figure 5. The solution on the axis at t=0.4 in problem 1 in sec- 
tion [47T] The 3 different reconstruction methods (minmod: blue, 
3-minmod: sky blue and PPM: magenta) are shown. 



suits minmod, it cannot produce the sharp features of the 
shock seen in the ID test: this is due to the dissipation 
caused by the imbalance between the net flux and the source 
term. We can also see that 3minmod and PPM yield quite 
similar results. We checked the differences in solutions at dif- 
ferent z = const, planes and found that they are negligibly 
small (~ 1CP 13 ) compared to the truncation errors. Over- 
all, however, although the 2D results show more dissipation 
than the ID ones, the relative differences in the solutions are 
not significant (for our purposes). In particular, both agree 
acceptably with the analytic forms. 



5 STATIONARY STAR TEST 

The tests just reported did not involve the effects of the 
gravitational field. In this section, we test our treatment of 
the Poisson equation for the gravitational potential as well 
as the hydrodynamics. 

With an ideal code, the evolution of a stationary star 
should also be stationary. However, in practice, all codes 
that dynamically evolve stationary states show some level 
of fluctuation due to finite grid resolution and intrinsic er- 
rors in the numerical scheme used. In this section, we show 
the time evolution of the physical quantities of non-rotating 
and the rotating stars, and investigate the dependence of 
this time behaviour by changing the resolution of the sim- 
ulations. Specifically, we use 3 different grid resolutions : 
65 x 65, 129 x 129 and 257 x 257, where 1/2 of the grid 
points span the star at the equator. 

Our initial models of rotating stars a re generated using 
Hachisu's Self-Consistent Field (HSCF: lHachisul Il986allbl) 
metho d — details of the procedure are described in lKim et alj 
(2009). In order to generate equilibrium models, we choose 
1) the maximum rest mass density, po 13 *) 2) the rotation pa- 
rameter, A, which describes the differential rotation and 3) 
the axis ratio which determines how fast the star is rotating. 
We must also specify the equation of state (EOS) in our con- 
struction of the initial model. Here, we used the polytropic 
EOS Eq. (HU) with K = 100 and TV = 1. We choose a maxi- 
mum density value of p^ ax = 1.28 x 10~ 3 , which, with this 
EOS, produces a 1.4M0Star in the non-rotating case. For 
the rotating models, we only consider rigid body rotation, 
which is obtained when we chose a very large value of A. 
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The axis ratio is specified to be 0.75 resulting in an orbital 
frequency of 6f 1Hz. 

Even with our use of the multigrid technique — which is 
generally an efficient method for solving elliptic equations — 
we still find solution of the Poisson equation for the gravi- 
tational potential to be computationally expensive. We thus 
calculate <3> only every 40 time steps to reduce the time spent 
in the Poisson solver, and find that this produces results 
which are nearly equivalent to those obtained when the Pois- 
son equation is solved at each time step. However, we use 
time-extrapolated values for the gravitational potential at 
the time steps between solves of the Poisson equation in or- 
der to avoid discontinuities in the primitive variables, when 
abrupt changes of the gravitational potential occur. We find 
that these discontinuities give rise to very unnatural dis- 
sipative effects in the simulation resulting, for example, in 
a rapid decay in the amplitude of radial oscillations, even 
when radial perturbations are explicitly introduced. 

Figure 6 shows the time evolution of the relative changes 
of the maximum density ([p? ax (t) - po lax (0)] /po**(0)) for 
non-rotating (left panel) and rigidly rotating (right panel) 
stars. For the stationary stars, we use the Cowling approx- 
imation, which assumes the gravitational potential is fixed. 
This gives efficient evolution of the stars, and can also be 
used as a testbed for fully coupled evolutions. The results 
computed using the Cowling approximation are shown in 
the top figures. The maximum density slowly increases with 
time for the rotating star while it decreases for the non- 
rotating star. For grid resolutions greater than 65 x 65 the 
rate of change is almost independent of resolution for the 
spherical star, but a slow decrease with resolution is seen 
for the rotating star. We define the following dimensionless 
rate of change: 



K 



dlnp n; 



dt 



(24) 



where we use id yn = 1 / v / Po lax f° r simplicity. We use 
this quantity — as computed from the highest resolution 
simulations — as a label in the figures. The values of 1Z are 
within 3 x f -7 for non-rotating star and are about 10 times 
larger for the rotating star, again with a maximum resolu- 
tion of 257 x 257. The inverse of 7Z can be interpreted as 
the time (in units of the dynamical time) that the simula- 
tion could be carried out until the results deviate from the 
true solution by 0(1). Our results indicate that the error 
would become ~ 1% in 30,000 and 3,000 dynamical times 
for non-rotating and rotating stars, respectively. We also 
carried out very long time simulations and found that 7Z 
becomes smaller even though it appears to be almost con- 
stant in the figures. From these results, we conclude that we 
can use the code to evolve stellar configurations for several 
thousand or more dynamical times. 

It is also very important to check the constancy of the 
conserved quantities with respect to simulation time. In our 
formulation we have 2 conserved quantities: the total rest 
mass, Mo, and the total angular momentum, J, which are 
computing using 



M = / DdV i3) = 2tt 



2tt 



PoW 



3/2 



(1 + 2$) 



RdRdZ, 



(1 + 2$) 



5/2 



R^dRdZ, 



(25) 
(26) 



-0.0001 



^ -0.0004 




Figure 7. The deviation of total rest mass (upper panel) and 
angular momentum (lower panel), which should remain constant, 
from their initial values with time, for the same models shown in 
Fig. 6. Results computed with three different resolutions (65 x 
65:red, 129 X 129:dark blue and 257 X 257:sky blue) are shown. 



respectively, and where dV^ denotes the 3-dimensional vol- 
ume element. Figure 7 shows the time evolution of these two 
conserved quantities: total rest mass (upper panel) and total 
angular momentum (lower panel). We show the results only 
from the rotating star since there is, of course, no angular 
momentum for non-rotating stars. The deviation of the total 
rest mass from the initial value has two features: short-term 
fluctuations and long term average behavior. The shot-term 
fluctuations depend on the grid resolution, but the average 
slopes are almost independent of the resolution. We label the 
graphs with 1Z m and TZj in a manner analogous to Eq. (|24|l 
and Fig. [5] and use these quantities to measure the long-term 
stability of the code. Their measured values are consistent 
with the ones for the central density (1Z). The behaviour of 
TZm is quite similar for the 3 different grid resolutions, but 
Rj shows considerable dependence on the grid resolution. 
We have seen above (see Figure 6) that the central den- 
sity fluctuation is significantly dependent on grid resolution 
only for the rotating models. We conclude that the main 
reason for this resolution-sensitive behavior is the fact that 
angular momentum conservation is sensitive to grid resolu- 
tion. Therefore, simulations for rotating stars require high 
grid resolution, otherwise angular momentum conservation 
will fail, and other stationary properties of the start (such as 
central density) will also show substantial, and non-physical, 
time evolution. 
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Figure 6. The time evolution of the maximum rest mass density changes( [pQ Ilax (t) — P™ ax (t = 0)] / p Ilax (t = 0)) with different resolutions 
(65 X 65: red, f29 X f29: dark blue and 257 X 257: sky blue). Top figure shows results when we fix the metric (Cowling approximation) 
while the bottom one shows the case where we consider the fully coupled dynamics. In the left panel, we show the figures for a spherical 
(non-rotating) star while the right panel shows the corresponding figures for a rigidly rotating star with axis ratio= 0.75, which give a 
rotational frequency of 6f fHz 



6 RADIAL PULSATION FREQUENCY TEST 

Even without any explicitly-added perturbations, it is nat- 
ural for our numerical simulation of stationary stars to give 
rise to normal mode oscillations due to intrinsic numerical 
errors. These errors occur for a variety of reasons, includ- 
ing 1) truncation error due to the discretization scheme, 
2) the artificial atmosphere (floor) whereby the primitive 
variables (pressure, density) are restricted from falling be- 
low minimum values to avoid code crashes (the sound ve- 
locity becomes unbounded when vacuum is encountered in 
the numerical calculations) , and 3) the numerical limitation 
in describing the stellar surface. Furthermore, the artificial 
atmosphere is known to excite higher overtone modes. 

The frequencies of various modes depend only on the 
structure of a given star, and can be calculated by various 
methods. As explained above, our stationary models oscil- 
late even when we do not explicitly introduce external or in- 
ternal perturbations. We attempted to compare the frequen- 
cies of the modes excited in our models with those obtained 
by normal mode analysis. The fundamental mode (F-mode 
hereafter) frequency is very closely related to the dynamical 
time (~ \/yfp) and the associated overtones have frequen- 
cies of similar order. 

Although using calculations based on cylindrical coor- 
dinates is not an efficient way to compute radial pulsations, 
our code should still be able to approximately compute the 
correct pulsation frequencies. The detailed perturbation for- 
mulations and numerical methods we use for investigating 



the radial pulsations are described in Appendix [X] For ini- 
tial conditions we use a non-rotating equilibrium star with a 
baryon mass 1.4Mq . we performed the test with and without 
the Cowling approximation, and In order to obtain the mode 
frequency from the simulations, we analyzed the fluctuation 
of the maximum density with time. 

Specifically, we carried out Fourier transformation 
on the maximum de nsity using the FFTW package 
l|Frigo fc Johnson! 120051 ). To obtain better resolution in the 
frequency domain, we use the zero-padding method which 
adds additional zeros at the end of the time series data, effec- 
tively using interpolation between points following the basic 
Fourier transformations. During the process of obtaining a 
frequency having a maximum sinusoidal amplitude, leakage 
may also cause additional errors. To reduce the effects of this 
leakage, we multiply the time series by a window function. 
Here we used the Hamming window function defined by 



• ; - I).". I + 0. Kims ( ^ J 



(27) 



where j is the index of the grid points and N is the total 
number of points, prior to zero-padding l|Harrislll978h . 

Although, as described above, some modes are excited 
simply due to numerical error, their amplitudes are too small 
to be accurately extracted from the simulation. We therefore 
introduce an explicit perturbation which can more strongly 
excite the radial modes. The perturbation that we used is 



Spo = B s sin ( 7r — 



(28) 
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Table 3. Comparison of mode frequencies obtained by numerical 
simulation (2nf) and by linear analysis (cr) 
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Figure 8. The mode amplitudes of maximum density as a func- 
tion of frequency of the star with a baryon mass 1.4M0. The 
vertical red dotted lines show the frequency of the radial pulsa- 
tion modes computed using the perturbation method. The top 
panel shows the result when we use the Cowling approximation, 
where the gravitational potential is assumed to be fixed. In the 
bottom figure, we obtain the gravitational potential every few 
time steps. In the both panels the 3 curves show results obtained 
using 3 different grid resolutions (sky blue: 257 X 257, dark blue: 
129 X 129, and red: 65 X 65) 



where B 3 is the perturbation amplitude which we set to 
B s = 0.001. 

Figure 8 shows the result after Fourier transformation 
of the time series data given by the differences in maximum 
density relative to the initial time (pcF^M - Po lax (i = 0)), 
and using calculations at different resolutions. For compar- 
ison purposes, the vertical red lines show the results com- 
puted from linear analysis. The mode labeled as F is the fun- 
damental mode, while H n denote the n-th overtone radial 
modes. The results shown in the figure can be summarized 
as follows: 

(i) The most excited mode with the perturbation given 
by Eq. (|28|) is the F-mode. By changing the nature of the 
perturbation we could make one of the overtones the most 
highly excited. 

(ii) At low resolution, the code cannot identify high fre- 
quency modes. The reason for this is the lack of spatial, not 
temporal, resolution. The eigenfunctions describing higher 
overtones have large gradients near the surface which cannot 
be accurately represented in the low resolution calculations. 

(iii) The frequency increases when we use the Cowling 
approximation. This is a well-known phenomena irrespective 
of whether Newtonian or general relativistic gravitation is 
used. This issue is discussed in more detail in the appendix. 



Table [6] shows the mode frequencies computed from lin- 
ear analysis as well as the numerical simulations. Again, the 
stellar model is a non-rotating spherical star of mass 1.4 
M©. The relative difference between the linear and full nu- 
merical results is listed in the last row. Here the numerical 
simulations have been carried out using the highest reso- 
lution (257 x 257), and we list results computed with and 
without the Cowling approximation. The frequencies we ob- 
tained from the numerical simulation with 257 x 257 grid res- 
olution have relative differences from those computed from 
linear analysis of at most 0.1%. We thus conclude that the 
radial mode frequencies computed from our code agree very 
well with the ones calculated from linear theory. The largest 
difference of 0.1% was found in the second overtone (mode 
H2), while for other modes we did not find any measurable 
difference. 



7 SUMMARY AND DISCUSSION 

We have developed a new hydrodynamical code which 
adopts a pseudo-Newtonian treatment of the gravitational 
field. This code uses the so called "Valencia formulation" for 
the hydrodynamical equations. From the computational per- 
spective, the code is modular and includes many reconstruc- 
tion schemes such as slope limiting techniques (minmod, 
MC, 3rd order minmod, etc.), PPM and ENO (WENO). 
In ID shock tube tests, we assessed code accuracy relative 
to analytic solutions and computed convergence rates of the 
errors. We found that the minmod method gives the most 
diffusive results, smoothing out complex features near dis- 
continuities. As a result it cannot be used to accurately de- 
scribe stellar surfaces, which are characterized by stiff den- 
sity changes. The MC method gives the most promising re- 
sult in the shock tube test and has second order accuracy. It 
can capture discontinuities very well in the pulsation mode 
test, but also yields additional non-physical effects such as 
the excitation of the higher order overtones near the stel- 
lar boundary. The 3minmod and PPM methods can provide 
higher order accuracy and we have found that they can also 
describe the stellar surface well. 

In the code we also implemented 3 different flux ap- 
proximation schemes: Roe, Marquina, and HLLE. Although 
the results in this paper were all computed using the HLLE 
approach — which is the most dissipative of the three — we 
have also found that for the simulations we have considered 
all produce very similar results. 

In the multigrid module for computing the gravitational 
potential we have implemented both second and fourth or- 
der finite-difference discretizations. The actual value of the 
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gravitational potential is slightly different if we change the 
order of accuracy. However, the changes of maximum density 
in time show very little sensitivity to the order of approx- 
imation, and we consider the difference between the use of 
the second or fourth order method to be insignificant. 

In the stationary star test which is described in sec- 
tion [5] we evolve equilibrium solutions describing both non- 
rotating and rotating stars using our code. Our code shows 
stable long-time behavior of the maximum density and con- 
served quantities. Based on the rates of change in the max- 
imum density, total mass and total angular momentum, we 
estimate that our code can be used to study evolution in 
excess of 3,000 dynamical times with 1% error. 

In the radial mode test described in section|6] modes are 
obtained from the Fourier transformation of the maximum 
density fluctuations. We also computed normal modes by 
linear analysis (see Appendix [A} and found that the mode 
frequencies generated by our code agree with the results 
from linear analysis almost perfectly (less than 0.14%). 

This code can be applied to the following astrophysical 
scenarios: 

(i) Phenomena associated with isolated rotating neutron 
stars, such as axisymmetric pulsations. Since our approach 
can be applied to mildly compact stars, it is very useful to 
determine the amplitudes and frequencies of the radial and 
non-radial modes. 

(ii) Accretion disks around a neutron star or black hole. It 
is not sufficient to treat a disk around a compact object using 
Newtonian gravity, since the gravitational field is not weak 
there. In addition, because the rotational velocity of the disk 
is a significant fraction of c, we should also take into account 
special relativity in our treatment of the hydrodynamics. 
Our code can be a very good tool for accretion disk studies. 

This work was supported by the NRF grant 2006-341- 
C00018, by NSERC, and by the CIFAR Cosmology and 
Gravity Program. 
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APPENDIX A: PERTURBATION EQUATION 

The eigenfrequencies and eigenfunctions of the radial pul- 
sation of stars are well-known in Newtonian hydrodynamics 
as well as in the general relativistic case. However, the cor- 
responding formulation has not been previously presented 
for our pseudo-Newtonian approach. Here, we describe the 
linearized equations that can be used to obtain eigenfre- 
quencies and eigenfunctions of the normal modes of spher- 
ical stars using this approximation, and following the gen- 
eral relativistic framework described in iMisner et al. 1973 
(MTW hereafter). First, to describe stellar oscillations — 
such as those occurring on the surface — it is much more 
practical to use a Lagrangian description rather than the 
Eulerian one adopted in section [2] The relation between the 
the Eulerian and Lagrangian perturbation is (see, e.g., Cox 
1980), 

Af(t,r) = 6f + tiC, (Al) 
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where £ is a Lagrangian variation in space. The law of baryon 
number conservation(V M (rm M ) = 0) gives 



An = 



-no 



> 2 ao(r 2 a 3 C)'-3a 1 <5a], 



(A2) 

where a — yl + 2$, n is the baryon number density 
and ' denotes differentiation with respect to r (See MTW 
Eq.(26.7)). The relation between n in Eq. (|A2|) and p is 
po = mi,n, where mt is baryon mass and the subscript 
denotes the unperturbed state. 

Another perturbation equation comes from the adia- 
batic equation of state which offers a much easier way to 
find the pressure variation: 

Since the Lagrangian variations commute with total differ- 
entiation (denoted by d), Eq. <|A3[) becomes 



r _ n AP 
P An ' 



(A4) 



In addition, Eqs. (|A"Tj) . (|A"2"1) and (|A"4)) give the following 
pressure variation equation: 

SP = -FP [r- 2 a 3 (r 2 ao 3 0' - 3a : 1 Sa} - C^o- (A5) 
The energy conservation equation (u,,V v r'"') gives 

(A6) 



Ap= Po + Po An 



n 

Note that po is the energy density in the unperturbed state, 
rather than the rest mass density used in the main text. 
Combining this with Eq. (|A2|) . we obtain the equation for 
the energy density variation 

Sp = -(po + Po)[r~ 2 a 3 o(r 2 ao 3 0' - 3a^5a] - Cp'o- (A7) 

The main difference here relative to the general rela- 
tivistic case arises in the computation of the perturbation of 
the gravitational potential. The Poisson equation gives 



— (aoSa)' + (aoSa)' 
r 



4tv(8p + 35P). 



(A8) 



Note that we should use only the Eulerian variation in this 
equation since Lagrangian variation does not commute with 
partial differentiation. Eq. (26.16) in MTW involves only 
first order differential equations — i.e. the second order dif- 
ferentiations are rewritten in terms of first order ones. On 
the other hand, in our case we cannot find any equations 
which can be used to eliminate the second order differenti- 
ation. That means that we need to find one more boundary 
condition to solve this equation. 

Finally, the equation of motion of the fluid is obtained 
from the 4-acceleration (a M = u"V„!i^), 

(po + Po)ao *C = -5P' - (5p + SP^a'o 

- (po + Po)(ao 1 8a - «o 2 a' Q 8a). (A9) 

Under the assumption of the adiabatic nature of the 
oscillation, normal modes are standing waves, and thus space 
and time variables can be separated as follows: 

<(r,t)=eMe lCTt . (AI0) 
Then, we can rewrite the equations using £ and a, 

8P = -TPo[r- 2 a 3 (r 2 ao 3 0' ~ Sa^Sa] - fPo" (All) 
8p = -(po + Po)[r- 2 a 3 o(r 2 ao 3 0' ~ Sa^Sa] - £p(AI2) 
(po + PoK" 4 <A = SP' + (5p + SP^a'o 

+(po + Po)(Q!q da - Qq 2 a 5oO\I3) 



To solve Eqs. (fATT) - (fAf3| . we need to impose appropriate 
boundary conditions. The first condition is that £/r should 
be regular at the origin, and the second one is that the pres- 
sure variation at the surface must vanish, i.e., 



- = finite at r — 0, 
r 

AP(r = r„) = 0. 



(A14) 
(A15) 



Unlike the general relativistic case, we cannot substitute Sa 
and 5a' in terms of other variations such as 5p and SP. 
Therefore, we need an additional boundary condition for Eq. 
(|A8|) . We use the properties of the gravitational potential to 
obtain extra conditions. First, from the condition that the 
gravitational potential should be regular at the center we 
obtain 



Sa' = at r — 0, 



(A16) 



Second because the gravitational potential should fall off as 
1/r beyond the stellar surface, we have 



5$' + £i=0. 
r 



(A17) 



When we apply the above equation at the stellar boundary 
(r = r s ), we get 



Sa' 



Sa 2 



2rSa 



at r 



(A18) 



Since Eqs. (|A11|) - (|A13|) and (|A8|) are coupled, we use 
an iterative method to solve them. 

For the case of the Cowling approximation, which as- 
sumes that the gravitational potential is fixed (Sa = 0), the 
equations simplify considerably: 

SP = -rPo[r- 2 a 3 (r 2 ao 3 0'] - tH 



Sp = -(po + Po)[r- 2 a 3 (r 2 ao 3 0'] - £Po 



(A19) 
(A20) 

(po + Po)a " V< = SP' + (Sp + 5P)ao 1 a' (A21) 

If we compare the above equations with Eqs. (|AI I[) - (|AI 3[) , 
we observe that every coefficient of Sa is negative: therefore, 
as mentioned in the main text, a increases when we apply 
the Cowling approximation. 

We show the solution for f/r for the 1.4Mq star with 
K — 100 and N = 1 with and without the Cowling approx- 
imation in Figure IA1I The a values corresponding to each 
mode are summarized in Table [6] which appears in the main 
text. 
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Figure Al. Radial pulsation cigcnfunction of a 1.4M0Star. The equation of state that we use is the polytropic one with K = 100 and 
N = 1. In this figure, § = f/r and £ s = £(r = r s ) where r s is the surface radius. The dashed and solid lines represent the results with 
and without the Cowling approximation, respectively. Each panel shows different modes (top- left (F), top-right (Hi), middle- left (H2), 
middle- right (Hz), bottom-left (H4) and bottom-right (Hr,)) which have different oscillation frequencies. 



